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We calculate the NRQCD matrix elements for the decays of the lowest-lying S- and P- 
wave states of charmonium and bottomonium in quenched lattice QCD. We also compute 
the one-loop relations between the lattice and continuum matrix elements. 

1. Heavy-Quarkonium Formalism 

Heavy-quarkonium systems are nonrelativistic. In the CM frame, the average 
quark velocity v satisfies v 2 <C 1, where v 2 0.3 for charmonium, and v 2 w 
0.1 for bottomonium. This fact allows one to describe heavy-quarkonium systems 
conveniently in terms of the effective field theory Nonrelativistic QCD (NRQCD) u 
NRQCD accurately describes processes in which pq < Mq. Its utility stems from 
the fact that it can be used to decouple short-distance (~ 1/Mq) processes from 
long-distance (~ quarkonium size ~ 1/(Mqv)) processes. 

QQ annihilation occurs at a distance of order 1/Mq, so NRQCD does not de- 
scribe the details of that process. In NRQCD, the short-distance part of the am- 
plitude for QQ — > light hadrons — > QQ is pointlike, and the entire amplitude is 
described by a four-fermion interaction. The quarkonium total annihilation rate is 
proportional to the imaginary part of the matrix element in the quarkonium state 
of the appropriate four-fermion operator. 

Coefficients of the four-fermion operators are determined by matching matrix 
elements in NRQCD to those in full QCD. The coefficients are short-distance quan- 
tities and, hence, are calculable in QCD perturbation theory. They are proportional 
to the IR-finitc parts of the parton-level annihilation rates. 

1.1. Factorization theorems 

Using these ideas, Bodwin, Braaten, and Lepage! have shown that a quarkonium 
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decay rate can be written as a sum of terms. Each term is the product of a long- 
distance matrix element of a four-fermion operator in the quarkonium state with a 
short-distance coefficient. For example, decay rates for S-wave quarkonia through 
next-to-leading order in v 2 are given by 

T( 2s+1 S,, -> X)=g 1 ( 2s+1 Sj) 2Imf 1 ( 2s+1 Sj)/M 2 

+F 1 ( 2s+1 Sj)2Img 1 ( 2s+1 S J )/M*. (1) 

Decay rates for P-wave quarkonia, to the lowest non-trivial order in v 2 , are given 

by 

Tf+ipj -> X)=H 1 ( 2s+1 Pj)2ImA( 2s+1 P J )/M* 

+Hs( 2s+1 Pj)2lmf s ( 2s+1 Sj)/M 2 . (2) 

The /'s and g's are the short-distance coefficients and are proportional to the rates 
for the annihilation of a QQ pair from the 2s+1 Lj state. Qx, Txi TC%, and TCs are 
the long-distance matrix elements, which we calculate in this paper. The subscripts 
1 and 8 indicate the color state of the QQ pair. (An octet state, in lowest order in 
v, consists of a QQg state.) 

1.2. Matrix elements 

The long-distance matrix elements are defined by 

^=(^1^/2) 6 X-X f (*/2) D i>\ l Px), 
H 8 =( 1 Pi|V t T°xx t T>| 1 P 1 ). (3) 

Here, tp and x^ are two-component spinors that annihilate a heavy quark and a 
heavy antiquark, respectively. 

The Qx and Hi terms in the decay rates (||) and (||) appear in the conven- 
tional, color-singlet modelB For the color-singlet matrix elements we can, to good 
approximation, take only the vacuum in the intermediate state.i Then we have 

Gi=\( 1 s \tfx\o)\ 2 (i + 0(« 4 )), 

Ty= Cs \^ X \0)(0\xHy Bf^So) (1 + 0(v*)), 
H l =\{ 1 Pi\^^ Dx|0)| 2 (l + O(« 4 )). (4) 
It then follows that, in the vacuum-saturation approximation, 

Si* f |i?5(0)| 2 , 
Z7T 

Hi« ^-\R' P (0)\\ (5) 
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where R(Q) is the radial wavefunction at the origin, and R'(0) is the derivative of 
the radial wavefunction at the origin. The matrix-elements (^) define a regularized 
i?(0) and a regularized R'(0) in QCD. The ratio T\jQ\ measures the average of p 2 
in the quarkonium state. 

The term in the P-wave decay rate (|^) that is proportional to 7is is absent in 
the color-singlet model. The matrix element Tig is proportional to the probability 
to find a QQg component in P-wave quarkonium, with the QQ pair in a relative 
S-wave, color-octet state. It is perhaps the most interesting of the matrix elements 
that we measure, since it corresponds to a field-theoretic effect of QCD and is, 
therefore, inaccessible through any potential model of quarkonium. 

2. Lattice Measurement of the Matrix Elements 

2.1. Euclidean- space formalism 

In Euclidean space, quarkonium two-point functions decay exponentially: 

lim (M(T)M + (0)) = lim {Q\Me~ HT M^\Q) = (0\M\l)e- E ' T (l\M^\0). (6) 

Here, M and M are quarkonium sources, and |^) is the lowest-lying state with the 
quarkonium quantum numbers. We deduce from (^|) that, in the vacuum-saturation 
approximation, the quarkonium matrix elements are proportional to the coefficients 
of e~ E ' T in the appropriate two-point functions. 

Similarly, a three-point function involving the four-fermion operator O has the 
behavior 

lim (M(T + T') 0{T) {0)) = (0\M\l) e - ElT ' (l\0\l}e~ E ' T (l\M^\0) 

T,T' — >oo 

=(M(T + T') Aft (0))(l\O\l). (7) 

Therefore, one can determine the expectation value of O in the quarkonium state 
by measuring the ratio of a three-point function to a two-point function. 
In practice, we rewrite the two-point function in (^) as follows: 

lim (M(T + T>) Aft(O)) - ^{M(T + T') M p (T)) (Mp(T) Aft(O)), (8) 

where Mp is a point-source interpolating operator, and C is the coefficient of e~ ElT 
in (Mp(T) Mp(0)). The sources in the three-point function match those in the two- 
point functions on the right side of (^). Hence, we can reduce noise by measuring 
the ratio of the three-point function to those two-point functions. The expectation 
value (i|C|Z) is then represented by the diagram in Fig. 1. 

2.2. Computational method 

We measured operator expectation values using noisy-point and noisy-gaussian 
sources and generating retarded and advanced quark propagators from each time 
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Fig. 1. Diagrammtic representation of (Z|0|Z). The large disks represent quarkonium sources (and 
sinks); the small disks represent quarkonium point sources and four-fermion operators. 

slice. We chose the Coulomb gauge for the field configurations. This gauge choice 
made the implementation of extended sources simpler, and it allowed us to replace 
covariant derivatives with normal derivatives, with errors of relative order v 2 . For 
matrix elements of covariant operators, we checked some of our results on non- 
gauge-fixed field configurations. 

In calculating heavy-quark propagators G(x, t) on the lattice, we used the non= 
relativistic formulation of Lepage, L. Magnea, Nakhleh, U. Magnea, HornbostelH 
We chose an evolution equation that is valid to the lowest non-trivial order in v 2 , 
that is, a lattice version of the inhomogeneous Schrodinger equation: 

G(x, t + 1) = (1 - H /2n) n ul t (l - ff /2n)"G(x, t) + d^ Q S t+h o, (9) 

with the initial condition G(x,i) = for t < 0. Here, H = -V (2) /(2M ) - ho, 
V^ 2 ) is the gauge-covariant discrete Laplacian, Mq is the bare heavy-quark mass, 
u o = ((l/3)Trt/pi aqU ette} 3 is the mean-field value of a gauge-field link, and ho — 
3(1 — uq)/Mq is the mean-field energy shift. This choice of ho is equivalent at 
leading order in v to tadpole improvement of the actional We chose n = 2. For our 
choices of the bare masses M , this is the minimum value of n for which the QQ 
propagators are free from lattice-artifact singularities. 

Note that J^i/Mq is suppressed by 0(v 2 ) relative to Q\. That is, this ratio is of 
the same order as terms that we have neglected in the evolution equation. However, 
in 3 Si — > Light Hadrons, 3 5i — > 7 + Light Hadrons, and 3 Si — > 37, the coefficient 
of T\/Mq is approximately —5 times that of ^iBB Thus, we feel that there is some 
merit in calculating J-'x/Mq, even in the presence of order v 2 errors. 

2.3. Parameters of the lattice simulation 

In our bottomonium and charmonium measurements, we used of 158 quenched 
gauge-field configurations on a 16 3 x32 lattice at 6/g 2 — 5.7. In the case of bottomo- 
nium, we also used 149 configurations on a 16 3 x 32 lattice at 6/g 2 — 6.0. We took 
our values of the bare heavy-quark masses from the determination by the NRQCD 
collaboration! M b0 = 1.5 at 6/g 2 = 6.0, M b0 = 2.7 at 6/g 2 = 5.7, and M c0 = 0.69 
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at 6/g 2 = 5.7. (Note that our definition of Mq is uq times that of the NRQCD 
collaboration.) We also used u = 0.87778701 at 6/g 2 = 6.0 and u = 0.8608261760 
at 6/g 2 = 5.7. 

In converting from lattice to physical units, we used the lattice spacings de- 
termined by the NRQCD collaboration. Their values are a -1 = 2.4 GeV for bot- 
tomonium at 6/g 2 — 6.0, a -1 = 1.37 GeV for bottomonium at 6/g 2 = 5.7, and 
a -1 = 1.23 GeV for charmonium at 6/g 2 — 5.7. At a given coupling, the lattice 
spacings are different for charmonium and bottomonium because the quenched ap- 
proximation leads to slight inconsistencies when one tries to use the physical spectra 
to fix the lattice spacing. 

2.4. Lattice results 

Our measurements revealed that the vacuum-saturation approximation is even 
more accurate than one would expect. For bottomonium at 6/g 2 = 6.0, corrections 
to the vacuum-saturation approximation for Gi are of relative size 1.3(1) x 10~ 3 . 
For charmonium at 6/g 2 = 5.7, the corrections are about 1%. The corrections to 
the vacuum-saturation approximation for H\ are also small. We assumed that the 
vacuum-saturation approximation is accurate for T\ as well. The numerical results 
that we present for Gi, F\, and Hi are the vacuum-saturation values. 

The results of our lattice measurements of the matrix elements are shown in 
Table 1. (The subscript L denotes lattice regularization.) The first error in each 



Tabic 1. Lattice values of the NRQCD matrix elements. 





charmonium 


bottomonium 


6/9 2 


5.7 


5.7 


6.0 


SlL 


0.1317(2) (12) 


0.9156(9)(65) 


0.1489(5)(12) 


J r iL(non)/g 1L 


1.2543(7) 


2.7456(8) 


1.3135(8) 


Fil{cov)/G\ L 


0.5950(5) 


2.1547(7) 


0.8522(5) 


3 r iL(non 2 )/GiL 


0.7534(4) 


1.2205(2) 


0.7775(5) 


Fil(cov 2 )/Gil 


0.5201(3) 


1.1111(2) 


0.6659(3) 




0.0208(2)(20) 




0.0145(6)(20) 




0.034(2)(8) 




0.0152(3)(20) 



quantity is statistical. Where two errors are given, the second error is an estimate 
of systematics associated with the parametrization of the functions used to fit to 
the propagators and with contamination from higher states. The arguments of Til 
indicate different lattice representations of the operator. The argument cov denotes 
a tadpole-improved naive gauge-covariant operator; the argument non denotes the 
simple, gauge-noncovariant, finite-difference operator in the Coulomb gauge. The 
subscript 2 indicates a difference operator with spacings of two lattice units. This 
softer lattice Laplacian was useful in controlling regulator-artifact power divergences 
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in the operator matrix element, which we shall discuss later. 

3. Relation of Lattice Matrix Elements to Continuum Matrix Elements 

At leading order in v 2 , the lattice and continuum matrix elements of the opera- 
tors that we measured are related as follows: 

0iL=(l + e)0i, 

^il = (1 +7)^1 + ^1, 

Wil = (1 + l)Hi + kHs, 

H SL =(l + v)tt 8 + (Hi. (10) 

(The continuum-regulated matrix elements have no subscript.) 

Note that, to leading order in v 2 , Gil has no T\ component. In fact, if one were 
to try to compute the addmixture of J-\, using the leading order NRQCD action 
that corresponds to (||), then the resulting inconsistencies in the treatment of v 2 
corrections would lead to uncanceled IR divergences. 

The coefficients in (|l(]) relate different UV regularizations of the operators. 
Therefore, they are short-distance quantities and, hence, are perturbatively cal- 
culable. The coefficients e, 7, 4>, t, 77 and ( are of order a s , while k is of order 
a 3 

We calculated these coefficients through order a s (one loop) in tadpole-improved 
perturbation theory.El First, we obtained analytic expressions for the integrands. By 
carrying out the time components of the loop integrations analytically, we were able 
to identify the IR divergent pieces, which are identical in the lattice and continuum 
matrix elements. After subtracting these divergent pieces, we evaluated the remain- 
ing IR-finite integrals numerically using VEGAS. 

Our numerical results (in lattice units) for the case of MS regularization of the 
continuum matrix elements are shown in Table 2. The accuracy of the coefficients of 
a s is better than 1%. The quantity £ depends at one-loop order on the factorization 
scale, which we took to be 1.3 GeV for charmonium and 4.3 GeV for bottomonium. 
These values are approximately equal to the MS c-quark and b-quark masses, 
respectively. 

Some of the coefficients of a s in the expressions for <f> appear to be large. How- 
ever, in physical units, Gi has dimensions (mass) 3 and T\ has dimensions (mass) 5 . 
Hence, we see from ( |l0| ) that (f> has dimensions of (mass) 2 . We can make <fi di- 
mensionless by dividing T\ and <f> by Mq. Similarly, we can make k and £ di- 
mensionless by dividing Hi and n by Mq and multiplying £ by Mq. Using the 
values, Mb = 5.0 GeV and M c = 1.5 GeV, we find that none of the dimensionless 
coefficients of a s is large. We conclude that the perturbation series is reasonably 
behaved in one-loop order. 

In setting the scale for a s , we made use of the method of Brodsky, Lepage 
and Mackenzie BIB For coefficients that arise from a positive integrand, the scale is 
about 1/a. In the case of integrands without definite sign, large cancellations can 
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Table 2. Coefficients relating lattice and continuum (MS) matrix elements. 





charmonium 


bottomonium 




5.7 


5.7 


6.0 


e 


-0.7326 a s 


0.2983 a s 


-0.4877 a s 


•y(non) 


-0.02578 a s 


-1.248 a s 


-0.9117 a s 


-y(cov) 


-2.860 a s 


-2.192 a s 


-2.560 a s 


■y(non 2 ) 


-0.2774 a s 


-1.096 a s 


-0.9236 a s 


<p(non) 


1.486 a 3 


10.90 a s 


4.418 a s 


<p(cov) 


0.3928 a s 


9.808 a s 


3.325 a s 


<p(non2) 


1.004 a 3 


6.096 a s 


2.863 a s 


L 


-0.7603 a s 


-1.852 a s 


-1.191 a s 


V 


0.09157 a s 


-0.03728 a s 


0.06096 a a 


c 


-0.1785 a s 


-0.006011 a a 


-0.01862 a s 



make the normalizing integral anomalously small and spoil the simplest scale-setting 
method. Therefore, we chose the scale to be 1/a for all of the coefficients, taking 
a s = a v (l/a) = 0.3552 at 6/g 2 = 5.7 and a s = a v (l/a) = 0.2467 at 6/g 2 = 6.0. 

4. Continuum Matrix Elements 

Substituting the lattice matrix elements and the lattice-to-continuum coefficients 
into ([To]) , we obtain the results for the continuum-regulated (MS) matrix elements 
shown in the first two columns of Table 3. The first two errors are the statistical and 
systematic errors from the lattice measurements. The third error is the systematic 
error from the neglect of terms of higher order in a s in the lattice-to-continuum 
coefficients. 

Errors from the omission of terms of higher order in v 2 in the evolution equation 
and in the operator mixing have not been reported in the first two columns of 
Table 3. We expect these errors to be of order 10% for bottomonium and 30% for 
charmonium. In the case of Gil, the NRQCD collaboration has given results that 
are accurate to next-to- leading order in v 2 . The weighted averages of the singlet- 
and triplet-state values are Gil = 0.133(4) for charmonium at 6/g 2 — 5.7 and 
Gil = 0.144(4) for bottomonium at 6/g 2 = 6.0, which are in good agreement with 
our results. This suggests that the approximate effect of the corrections of higher 
order in v 2 is to split the values of the matrix elements for a multiplet of spin states, 
without changing their spin average. 

The second column of Table 3 does not include errors that arise from uncertain- 
ties in the physical values of a -1 . We estimate, from the results of the NRQCD 
collaboration, that these errors are 7% for Gi in charmonium, 13% for Ti.\ in char- 
monium, 13% for Gi in bottomonium, and 23% for 7ii in bottomonium. 

In addition, there are errors associated with the quenched approximation, for 
which we have no quantitative estimate. 
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Table 3. Continuum-regulated (MS) matrix elements. 





lattice 


experiment 




lattice units 


physical units 




charmon 
Si 

Fi/Qi 
Hi 

Hs/Hi 


um 6/g 2 = 5.7 
0.1780(3)(16)(± 3 f 9 ) 

0.05 — 0.54 
0.0285(2)(27)(±42) 
0.086(l)(6)(lg) 


0.3312(6) (30) (t%l) GeV 3 

0.07 — 0.82 GeV 2 
0.0802(6) (77) (till) Gey5 
0.057(l)(4)(± 2 ]j GeV" 2 


0.36(3) GeV 3 
0.057 GeV 2 
0.077(19) (28) GeV 5 
0.095(31)(34) GeV -2 


bottomo 

Si 


lium 6/g 2 = 5.7 
0.8279(8)(59)(±i°f) 
-3.7 — 0.2 


2.129(2)(15)(± 2 ^) GeV 3 
-6.9 — 0.4 GeV 2 


3.55(8) GeV 3 


bottomo 

Qi 

F1/Q1 
Hi 

Hg/Hi 


lium 6/g 2 = 6.0 
0.1692(6)(14)(+^6) 

-0.34 — 0.28 
0.0205(9)(28)(1^) 
0.0151(2)(14)(± 33 ) 


2.340(8)(19)(+^ 3 ) GeV 3 
-2.0 — 1.6 GeV 2 
1.63(7) (23) (tig) GeV 5 
0.00262(3) (24) (l^)GeV" 2 


3.55(8) GeV 3 



4.1. Experimental values of the matrix elements 

The third column in Table 3 gives phenomenological results for the matrix ele- 
ments. Qi was extracted from the measured decay rates for J/tp — ► e + e~, r\ c — > 77 
and T -> e+e- (Ref. using the expressions in given in Eef. |. The value for 
F1/G1 for J/ip is from the calculation of Ko, Lee and Song.til Tti and Tis/Tti for 
Xc are from Ref. [l^. There is no published data for \b decays into light hadrons, 
photons, and/or leptons. To extract the phenomenological matrix elements for Q\ 
and Hi, we used the values M 6 (pole) = 5.0 GeV (Ref. g), M c (pole) =1.5 GeV 
(Ref. |l|), a s (M c ) = 0.243, a s {M b ) = 0.179, a(M c ) = 1/133.3, and a(M b ) = 1/132. 

In the third column of Table 3, the first error is experimental, and the second, 
where it is given, is theoretical. Where no theoretical error is given, it is at least 
as large as the uncertainty from the neglect of terms of higher order in a s in the 
calculation of the short-distance coefficients. This uncertainty is of nominal size 25% 
for charmonium and 20% for bottomonium. Errors that arises from uncertainties in 
the heavy-quark masses have not been reported in the third column. The NRQCD 
collaboration quotes an error of 4% for the 6-quark mass, which leads to errors of 
8% in Qi and 16% in Hi. In the case of the c-quark mass there is, as yet, no reliable 
determination from a lattice calculation. 

5. Discussion 

In the case of charmonium, our results for Qi, Hi, and H$/Hi are in agreement 
with experiment, but the errors are large. It is interesting to note that this agree- 
ment would have failed if we hadn't included the lattice-to-continuum corrections 
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to the matrix elements. 

The ratio T\jQ\ is poorly determined, largely because of uncertainties in the 
coefficient <fi that gives the mixing of T\ into Q\ . The mixing of T\ into Gi is power 
UV divergent. Since the mixing is UV dominated and begins at one loop, we expect 
it to be of order a s . On the other hand, in the continuum, J 7 i/(M 2 Q 1 ) is of order 
v 2 -c 1. Therefore the effect of mixing on the ratio T\j (M 2 Q\) is of relative order 
a s /v 2 and is large for both charmonium and bottomonium. 

Nevertheless, we can conclude that JFi/(M 2 Qi) is no larger than 0(v 2 ), in agree- 
ment with the NRQCD scaling rules.tM T\jQ\ is probably positive for charmonium 
and negative for bottomonium. Note that, because the continuum matrix element 
T is gotten by subtracting UV divergences, it need not be positive. 

For bottomonium, the lattice result for Q\ is 35 — 40% below the experimental 
value. We know from the results of the NRQCD collaboration that at least part of 
this discrepancy is due to the quenched approximation.Lrli-3 There is good agreement 
between our results at 6/g 2 — 5.7 and Q/g 2 = 6.0, which confirms the expected 
renormalization-group scaling behavior. 

Our results for the P-wave matrix elements for bottomonium can be translated 
immediately into predictions for bottomonium decay ratesEl The values for indi- 
vidual matrix are probably subject to large corrections from the quenched approxi- 
mation, but the ratio of octet to singlet matrix elements may be less susceptible to 
this source of error. 

Aside from quenching, the largest uncertainties in the matrix elements come 
from neglect higher-order (in a s ) corrections to the lattice-to-continuum coefficients. 
One might remedy this situation by using lattice methods to compute the relations 
between the lattice matrix elements and the momentum-subtracted continuum ma- 
trix elements nonperturbatively, as suggested by Martinelli and Sachrajda.ta The 
momentum-subtracted matrix elements could then be converted to MS matrix el- 
ements in continuum perturbation theory. 

In the continuum, ALS regularization of the operator matrix elements leads to 
renormalon ambiguities Jl3 These ambiguities are of the same order in v 2 as the ma- 
trix elements of operators of higher dimension. In the case of Hg, we expect such 
ambiguities to be small, since H\ first mixes with Hg in order a^. Renormalon am- 
biguities are absent in the case of hard-cutoff regulators, such as the lattice. That 
is, they are an artifact of the regulator (factorization) scheme that one chooses 
to define NRQCD. The consistency of NRQCD as an effective theory guarantees 
that regulator-scheme dependence is absent in physical quantities. Hence, renor- 
malon ambiguities cancel in decay rates if one computes the NRQCD short-distance 
coefficients and the lattice-to-continuum coefficients to the same order in a s . 

It is interesting that, for both charmonium and bottomonium, the values of 
Tis/TCi that we obtain are in agreement with a crude phenomenology^ In this 
phenomenology, one obtains Hs by solving the one-loop evolution equation for Wg, 
under the assumption that Hs vanishes below a scale Mqv. The one-loop evolution 
of the decay matrix element Tig is the same as for the corresponding production 
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matrix element Ti. s . This suggests that 7i' s ks H.%. The production matrix element 
Ti.' 8 can be extracted from CDF data for charmonium productionlij and from recent 
CLEO data.0 Using our value for Hi, we obtain WjUx = 0.042(19) GeV -2 and 
H' 8 /Hi = 0.046(28) GeV" 2 , respectively, both of which are in good agreement with 
our result for Hs/Hi- 
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